Data

gaza_geno <- read.csv("D:/R_Gaza/gaza/data/killed-in-gaza.csv")
class(gaza_geno)
## [1] "data.frame"

Libraries

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.4.2

First data exploration

gaza_geno <- as_tibble(gaza_geno)
class(gaza_geno)
## [1] "tbl_df"     "tbl"        "data.frame"
gaza_geno
## # A tibble: 34,344 × 7
##    id          name                       en_name         age dob   sex   source
##    <chr>       <chr>                      <chr>         <int> <chr> <chr> <chr> 
##  1 طفللللللل   بنت اسيا نزار محمد ابووردة Bint Asaia N…     0 2024… f     h     
##  2 طفللللل     ابن ولاء     الدوواسة      Abn Wlaa Ald…     0 2024… m     h     
##  3 طفل رضيع66  أمل محمد جمعة أبو صويص     Amal Muhamma…     0 2023… f     c     
##  4 طفل رضيع6   ياسر صلاح ياسر الدلو       Yasr Salah Y…     0 2023… m     h     
##  5 طفل رضيع55  كارمل علاء وليد حمدان      Karml Alaaa …     0 2023… f     c     
##  6 طفل رضيع5   ميساء ديب محمود السكافى    Maisaa Diab …     0 2023… f     c     
##  7 طفل رضيع4   محمد نعيم نصر عياد         Muhammad Nai…     0 2023… m     c     
##  8 طفل رضيع3   كنان شادي هاشم مشتهى       Knan Shadi H…     0 2023… m     c     
##  9 طفل رضيع226 ابن/ انوار احمد محمد درويش Abn/ Anwar A…     0 2023… m     h     
## 10 طفل رضيع225 بنت/زينب محمد العبد نواس   Daughter of …     0 2023… f     h     
## # ℹ 34,334 more rows
death_by_age <- 
  gaza_geno %>% 
  group_by(age = as.integer(age)) %>% 
  summarise(death = n())

peak_10_ages <- death_by_age %>% 
  top_n(n = 10, death) %>% 
  pull(age)

age_quantiles <- quantile(gaza_geno$age)
age_quantiles
##   0%  25%  50%  75% 100% 
##    0   14   26   39  101
peak_10_ages
##  [1]  0 16 17 23 28 29 30 31 32 33

Plots

theme_set(theme_bw())

gaza_geno %>% 
  ggplot(aes(x = age))+
  geom_histogram(fill = "#4D5052", colour = "black", bins = 20, binwidth = 2)+
  scale_x_continuous(breaks = seq(0, 100, 20))+
  labs(title = "Gaza genocide - Death by Age", x = "Age", y = "Frequency of Deaths", 
       caption = "The frequency of murdered palestinians by age (murdered by the Israelis)")+
  theme(plot.title = element_text(face = "bold"), 
        plot.caption = element_text(face = "italic", colour = "#4D5052"), 
        plot.caption.position = "plot")

ggsave(filename = "geno_hist.pdf", path = "D:/R_Gaza/pal_map/images")
## Saving 7 x 5 in image
death_by_age %>% 
  ggplot(aes(x = age, y = death))+
  geom_line()+
  geom_segment(x = 14, xend= 14, y = 0, yend = 659, linetype = "dashed")+
  geom_segment(x = 26, xend= 26, y = 0, yend = 652, linetype = "dashed")+
  geom_segment(x = 39, xend= 39, y = 0, yend = 469, linetype = "dashed")+
  scale_x_continuous(breaks = seq(0, 100, 20))+
  labs(title = "Gaza genocide - Death by Age", x = "Age", y = "Number of Deaths", 
       caption = "The number of murdered palestinians by age (murdered by the Israelis)")+
  theme(plot.title = element_text(face = "bold"), 
        plot.caption = element_text(face = "italic", colour = "#4D5052"), 
        plot.caption.position = "plot")

ggsave(filename = "geno_linechart.pdf", path = "D:/R_Gaza/pal_map/images")
## Saving 7 x 5 in image

Maps

Setting the data

# Get palestinian cities, convert it to sf dataframe, assign projection to the coordinates and save it
library(maps) # Get palestinian cities
## Warning: package 'maps' was built under R version 4.4.2
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
library(sf)   # convert it to sf dataframe, assign projection to the coordinates and save it
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
data(world.cities)
palestine_cities <- world.cities %>% filter(country.etc == "Palestine")
class(palestine_cities)
## [1] "data.frame"
palestine_cities <- as_tibble(palestine_cities)
write_excel_csv(palestine_cities, file = "D:/R_Gaza/pal_map/data/palestinians_cities.csv")

names(palestine_cities)    #get the names for the coordinate variables
## [1] "name"        "country.etc" "pop"         "lat"         "long"       
## [6] "capital"
palestine_cities_sf <- st_as_sf(palestine_cities, coords = c("lat", "long")) #convert to sf df tibble
class(palestine_cities_sf) #check if class has changed to sf dataframe tibble
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
st_crs(palestine_cities_sf) #check for geographic projections
## Coordinate Reference System: NA
st_crs(palestine_cities_sf) <- 4326  #assign geographic projection (4326, i.e., longlat WGS84)
st_crs(palestine_cities_sf) #check if the geographic projection was assigned
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
st_crs(palestine_cities_sf)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
palestine_cities_sf #check if "lat" and "long" variables were changed to "geometry" variable.
## Simple feature collection with 244 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 31.27 ymin: 34.25 xmax: 32.55 ymax: 35.45
## Geodetic CRS:  WGS 84
## # A tibble: 244 × 5
##    name               country.etc   pop capital      geometry
##  * <chr>              <chr>       <int>   <int>   <POINT [°]>
##  1 'Abasan al-Jadidah Palestine    5629       0 (31.31 34.34)
##  2 'Abasan al-Kabirah Palestine   18999       0 (31.32 34.35)
##  3 'Abud              Palestine    2456       0 (32.03 35.07)
##  4 'Abwein            Palestine    3434       0  (32.03 35.2)
##  5 'Ajjah             Palestine    5143       0  (32.37 35.2)
##  6 'Anabta            Palestine    7311       0 (32.31 35.12)
##  7 'Anata             Palestine    9149       0 (31.81 35.28)
##  8 'Anin              Palestine    3717       0  (32.5 35.17)
##  9 'Anzah             Palestine    2006       0 (32.37 35.22)
## 10 'Aqbah Jabbar      Palestine    6337       0 (31.83 35.45)
## # ℹ 234 more rows
st_write(palestine_cities_sf, "D:/R_Gaza/pal_map/data/palestinian_cities", 
         driver = "ESRI Shapefile", delete_layer = TRUE)  # save the sf data frame as shapefile
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `palestinian_cities' using driver `ESRI Shapefile'
## Writing layer `palestinian_cities' to data source 
##   `D:/R_Gaza/pal_map/data/palestinian_cities' using driver `ESRI Shapefile'
## Writing 244 features with 4 fields and geometry type Point.
unzip(zipfile =  "D:/R_Gaza/gaza/data/localities_palestine_sf.zip" , exdir ="data/localities_palestine")
localities_palestine_sf <- st_read("data/localities_palestine")
## Multiple layers are present in data source D:\R_Gaza\pal_map\data\localities_palestine, reading layer `%D8%AA%D8%AC%D9%85%D8%B9%D8%A7%D8%AA_%D9%81%D9%84%D8%B3%D8%B7%D9%8A%D9%86___Localities_Palestine'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.
## Reading layer `%D8%AA%D8%AC%D9%85%D8%B9%D8%A7%D8%AA_%D9%81%D9%84%D8%B3%D8%B7%D9%8A%D9%86___Localities_Palestine' from data source `D:\R_Gaza\pal_map\data\localities_palestine' using driver `ESRI Shapefile'
## Simple feature collection with 619 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 3809204 ymin: 3661924 xmax: 3960023 ymax: 3836009
## Projected CRS: WGS 84 / Pseudo-Mercator
localities_palestine_sf
## Simple feature collection with 619 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 3809204 ymin: 3661924 xmax: 3960023 ymax: 3836009
## Projected CRS: WGS 84 / Pseudo-Mercator
## First 10 features:
##    OBJECTID REGIONCODE DISTCODE GOVCODE  LOCCODE                    NAMEAR
## 1         1          2       24    2470 24703430             عبسان الجديدة
## 2         2          2       24    2470 24703470                     خزاعة
## 3         3          2       24    2465 24653200            مخيم دير البلح
## 4         4          2       24    2465 24653275               وادي السلقا
## 5         5          2       24    2455 24552681 أم النصر (القرية البدوية)
## 6         6          2       24    2465 24653065             مخيم النصيرات
## 7         7          2       24    2470 24703420                   خانيونس
## 8         8          2       24    2470 24703425                 بني سهيلا
## 9         9          2       24    2470 24703445             عبسان الكبيرة
## 10       10          2       24    2470 24703485                   الفخاري
##                                  NAMEEN
## 1                     'Abasan al Jadida
## 2                               Khuza'a
## 3                    Deir al Balah Camp
## 4                         Wadi as Salqa
## 5  Um Al-Nnaser (Al Qaraya al Badawiya)
## 6                      An Nuseirat Camp
## 7                            Khan Yunis
## 8                          Bani Suheila
## 9                     'Abasan al Kabira
## 10                         Al Fukhkhari
##                                                                                     NOTES
## 1                                                                                    <NA>
## 2                                                                                    <NA>
## 3                                                                                    <NA>
## 4                                                                                    <NA>
## 5                                                                                    <NA>
## 6                                                                                    <NA>
## 7  تشمل المَواصِي (خانيونس)، وقِيزَان النَجَّار، وقَاعْ الخَرَابَة، وقَاعْ القُرِين، وأُمُّ كَمِيل، وأُّمُّ الكِلاب
## 8                                                                                    <NA>
## 9                                                                                    <NA>
## 10                                                                                   <NA>
##    Shape__Are Shape__Len                       geometry
## 1   4919436.7  11010.191 MULTIPOLYGON (((3824072 367...
## 2   9193057.5  14360.481 MULTIPOLYGON (((3826180 367...
## 3    337604.1   3122.714 MULTIPOLYGON (((3823061 368...
## 4   8752155.1  12611.773 MULTIPOLYGON (((3827289 368...
## 5   4306257.9   9547.706 MULTIPOLYGON (((3844490 370...
## 6   1099596.7  10312.221 MULTIPOLYGON (((3827522 369...
## 7  72893548.1  53756.789 MULTIPOLYGON (((3821981 367...
## 8   9181044.3  16897.216 MULTIPOLYGON (((3823716 367...
## 9  17770337.9  20809.409 MULTIPOLYGON (((3823887 367...
## 10 12900017.9  17496.895 MULTIPOLYGON (((3820790 367...
class(localities_palestine_sf)
## [1] "sf"         "data.frame"
st_crs(localities_palestine_sf)$proj4string
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
localities_palestine_WGS84 <- st_transform(localities_palestine_sf, crs = 4326)
st_crs(localities_palestine_WGS84)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
st_write(localities_palestine_WGS84, "D:/R_Gaza/pal_map/data/localities_palestine", 
         driver = "ESRI Shapefile", delete_layer = TRUE)
## Deleting layer `localities_palestine' using driver `ESRI Shapefile'
## Writing layer `localities_palestine' to data source 
##   `D:/R_Gaza/pal_map/data/localities_palestine' using driver `ESRI Shapefile'
## Writing 619 features with 10 fields and geometry type Multi Polygon.
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: One or several characters
## couldn't be converted correctly from UTF-8 to ISO-8859-1.  This warning will
## not be emitted anymore.
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 147785534.120117009 of
## field Shape__Are of feature 57 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 263963660.09277299 of
## field Shape__Are of feature 83 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 135039999.595703006 of
## field Shape__Are of feature 232 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 249401566.34472701 of
## field Shape__Are of feature 323 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 248853481.261718988 of
## field Shape__Are of feature 366 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 300296232.731445014 of
## field Shape__Are of feature 495 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 140128653.509766012 of
## field Shape__Are of feature 497 not successfully written. Possibly due to too
## larger number with respect to field width
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: Value 100428703.091796994 of
## field Shape__Are of feature 548 not successfully written. Possibly due to too
## larger number with respect to field width
st_crs(localities_palestine_WGS84)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
st_crs(palestine_cities_sf)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"

Creating maps

ggplot(localities_palestine_WGS84)+
  geom_sf(colour = "black")

library(leaflet)
library(leaflet.providers)
## Warning: package 'leaflet.providers' was built under R version 4.4.2
library(htmlwidgets)

col_fun <- colorFactor(c("#337bd7", "#c8143b"), domain = NULL)
itamar <- st_point(c(34.82375, 31.89530))
itamarpopup <- c("Where I.J.R was born")   
pal_popup <- paste0("<strong>City Name:</strong>", palestine_cities$name)


palestine_map <- leaflet(localities_palestine_WGS84) %>% 
  addPolygons(stroke = TRUE, weight = 0.7, color = "black", 
              fillColor = ~col_fun(REGIONCODE), 
              group = "Regions", fillOpacity = 0.8) %>% 
  addCircleMarkers(lng = 34.82375 , lat = 31.89530, group = "My House",
                   popup = itamarpopup) %>% 
  addCircleMarkers(data = palestine_cities, group = "Cities", popup = pal_popup, 
                   radius = 0.01, color = "#2c120c") %>% 
  addTiles(group = "OSM") %>% 
  addProviderTiles(provider = "Stadia.OSMBright", group = "OSMBright") %>% 
  addProviderTiles(provider = "Stadia", group = "Stadia") %>%
  addProviderTiles(provider = "Stadia.Outdoors", group = "Outdoors") %>% 
  addProviderTiles(provider = "Stadia.StamenTonerLite", group = "TonerLite") %>% 
  addProviderTiles(provider = "CartoDB", group = "CartoDB") %>% 
  addProviderTiles(provider = "CartoDB.Voyager", group = "CartoDB.Voyager") %>% 
  addProviderTiles(provider = "Esri", group = "Esri") %>% 
  addLayersControl(baseGroups = c("OSM", "OSMBright", "Stadia", "Outdoors", "TonerLite", 
                                  "CartoDB", "CartoDB.Voyager", "Esri"), 
                   overlayGroups = c("Regions", "Where I Live", "Cities"))
## Assuming "long" and "lat" are longitude and latitude, respectively
palestine_map